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Using self-energy-corrected density functional theory (DFT) and a coherent scattering-state ap- 
proach, we explain current-voltage (IV) measurements of four pyridine-Au and amine-Au linked 
molecular junctions with quantitative accuracy. Parameter-free many-electron self-energy correc- 
tions to DFT Kohn-Sham eigenvalues are demonstrated to lead to excellent agreement with exper- 
iments at finite bias, improving upon order-of-magnitude errors in currents obtained with standard 
DFT approaches. We further propose an approximate route for prediction of quantitative IV char- 
acteristics for both symmetric and asymmetric molecular junctions based on linear response theory 
and knowledge of the Stark shifts of junction resonance energies. Our work demonstrates that a 
quantitative, computationally inexpensive description of coherent transport in molecular junctions 
is readily achievable, enabling new understanding and control of charge transport properties of 
molecular-scale interfaces at large bias voltages. 



There is significant interest in using organic molecules 
as active components in next-generation energy conver- 
sion devices such as organic photovoltaics [1 and dye- 
sensitized solar cells [2 , where a roadblock to higher 
efficiencies is quantitative understanding and control of 
charge transport phenomena at interfaces. Understand- 
ing electronic energy level alignment and transport at 
interfaces between active organic layers and conducting 
electrodes has been particularly challenging [3, 4 . How- 
ever, recent scanning tunneling microscope-based break- 
junction (STM-BJ) experiments [5l|6] of molecular junc- 
tions -devices formed by trapping organic molecules be- 
tween macroscopic metallic electrodes- have reported ro- 
bust conductance [SHS], thermopower [TOHIl], switching 
behavior [TS'-'TF, quantum interference effects [TBHIH], 
spin- filtering phenomena p!9H2T] . and even full nonlin- 
ear IV characteristics [161 |22j [23] , establishing such junc- 
tions as unique and revealing windows into the physics 
of charge transport at the molecular scale. 

Given the diversity of known synthesizable organic 
molecules, and a lack of intuition connecting transport 
phenomena to a specific molecule and interface, a rig- 
orous yet pragmatic quantitative theory capable of pre- 
dicting charge transport properties of molecular junc- 
tions with chemical specificity is needed. Charge trans- 
port calculations require solution of the electronic struc- 
ture for an open system out of equilibrium [24lj29j . The 
majority of prior theoretical studies have relied on a 
Landauer approach, simplified to treat electronic inter- 
actions at a mean-field level within density functional 
theory (DFT) using either Green's functions [3QH33] or 
scattering-states [34]. However, within common approxi- 
mations to DFT, such as generalized gradient approxima- 
tions (GGAs) and hybrid functionals, Kohn-Sham orbital 
energies are known to yield unsatisfactory level align- 
ment [35ti4Q] , leading to overestimated conductance [241 - 
129] and even incorrect trends [41 relative to experiment. 



Improved treatment of exchange and correlation ef- 
fects within the junction [29l [37l [38l I4TH5Q] has been 
shown to lead to better agreement with experiment. In 
particular, recent calculations based on many-body per- 
turbation theory within the GW approximation [TSl \42\ - 
l44l |5T1 have been demonstrated to amend junction level 
alignment and have resulted in conductance in quantita- 
tive agreement with measured values [15] [421444] . How- 
ever, these approaches rely on computationally expen- 
sive methods ^371 |38l l45ti5Q] , and their accuracy come at 
the cost of their tractability for junctions relevant to ex- 
periments; moreover, direct comparisons to experiment 
have been limited to low-bias measurements in the lin- 
ear response regime. The ability to measure, reliably, 
conductance and nonlinear IV characteristics, and the 
promise of exploring new phenomena at higher bias, calls 
for quantitative computational studies of truly nonequi- 
librium steady-states. 

In this Letter, we demonstrate and apply a quan- 
titative framework capable of explaining measured 
IV characteristics for four molecules -4,4' diamino- 
stilbene (DAS), bis- (4- aminophenyl) acetylene (APA), 
1,6-hexanediamine (HDA), and 4,4'-bipyridine (BP)- in 
contact with gold electrodes for biases as high as IV. 
Whereas the magnitude of the currents and shape of the 
measured IV characteristic are not reproduced by a stan- 
dard finite-bias DFT approach [3QH34] , a method based 
on the GW approximation [51], successful for low-bias 
conductance and extended here to finite bias, leads to 
excellent agreement with experiment. Our approach en- 
ables inexpensive and accurate calculations of coherent 
transport properties, which can be used for the design of 
functional molecular junctions with nonlinear IV charac- 
teristics. 

Previous studies [151 1421444] have shown amine- and 
pyridine-bonded molecules can bind preferentially to un- 
dercoordinated gold atoms. Building on those studies, we 



2 




Voltage [V] Voltage [V] Voltage [V] Voltage [V] 

FIG. 1: I{V) characteristics for (a) DAS, (b) APA, (c) HDA, and (d) BP junctions. Top: Experiment (color map), DFT-PBE 
(green line), and DFT-PBE(V) (dark-green points). Bottom: Experiment (color map), DFT+E (grey line), and DFT(V)+i; 
(black points) . Error bars are added to the computed currents to indicate the spread associated with the three different contact 
geometries used here (see SI). Measured IV characteristics for DAS, APA, and HDA junctions are adapted from [22]. Lewis 
structures of each molecule at top. 



construct three geometries for each junction with trimer 
(3 gold atoms), trimer-adatom, and adatom binding mo- 
tifs (see SI for geometries). We relax all junctions using 
DFT within the GGA of Perdew, Burke, and Ernzenhof 
(PBE) [52] and a double-^-basis set as implemented in 
SIESTA [51 |54]. Details of our DFT calculations are 
provided in previous work [42l-[44j. Atomic positions are 
relaxed until Hellmann-Feynman forces are smaller than 
0.04 eV/A. We model our system by 7-layers of 16 gold 
atoms on both sides, the last 4 layers being constrained 
to the bulk geometry. Initial trial geometries are adapted 
from previous works |T5l I42H441 155) l56] . We consider BP 
junctions in a "high conductance" configuration (HighG), 
following the notation of Ref. [15^. 

In what follows, the conductance and IV characteristics 
are computed from a Landauer-like formula using a co- 
herent scattering-state approach, and include exchange 
and correlation contributions that correct for zero-bias 
junction level alignment, as explained below. In this 
framework, the IV characteristic is expressed in terms 
of an energy- and bias-dependent transmission function 
T(co'; V) as 
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where / is the lead Fermi-Dirac occupation function and 
uj has units of energy. The transmission function is com- 
puted as T {uj; V) = tr [t {u; V) {u; V)] , where t and 
are the transmission coefficients of the junction scatter- 
ing states computed as a function of bias [34l |57]. The 
zero-bias conductance is determined via linear response, 
as S = §,\y^Q = ^r{EF;V = 0), where Ep is the 
junction Fermi level at zero bias. 

We obtain the self-consistent steady-state density ma- 
trix from DFT-PBE as described in ^341^3 using an 8 x 8 



/c//-mesh and in a manner equivalent to prior work [ 30] - 
34 . For finite bias calculations, referred to as DFT(V), 
the density matrix includes a real-axis integration of the 
scattering-states in the bias window on an adaptive en- 
ergy grid [34 with a resolution of up to 10~^eV in 
the vicinity of the molecular resonances. The transmis- 
sion function is generated in a subsequent step using a 
16 X 16/c//-mesh. 

To correct for inaccuracies associated with DFT-PBE 
Kohn-Sham eigenvalues for quasiparticle energy level 
alignment, we employ a physically- motivated electron 
self-energy correction to the molecular orbital energies in 
the junction, DFT+S, following Ref. [42H44] . Formally 
derivable in the weak coupling limit of the GW approxi- 
mation [44l|58l|59], this adjustable-parameter-free model 
self-energy acts on the molecular subspace and consists 
of two terms: i) a gas-phase correction accounting for the 
difference between DFT highest occupied molecular or- 
bital (HOMO) and lowest unoccupied molecular orbital 
(LUMO) energies and their gas-phase ionization poten- 
tial (IP) and electron affinity (EA) (Sqp); and ii) an "im- 
age charge" term, accounting for the polarization energy 
associated with static non-local correlations between the 
electrons (or holes) on the molecule and in the metal that 
close the gap of the molecule upon absorption (Scorr)- 
In this work, we calculate Hqf using a ASCF method, 
and following previous studies [111 I42H441 [5T] , we approx- 
imate ScoRR with an electrostatic image charge model 
(see SI for details). Since these self-energy corrections are 
large compared to those for the metallic bulk and surface 
Au states [60], especially for states near E^?, we neglect 
similar corrections to the Au states. Values for Hqp and 
ScoRR foi" each junction appear in the SI. This approach 
was shown to lead to quantitative agreement with pho- 
toemission experiments regarding the level alignment of 
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Junction Exp DFT-PBE DFT+E 



DAS 


10 


108 


- 125 


5.7-5.9 


APA 


8 


75 


-83 


6.0-6.4 


HDA 


1.2 


5.1 


-5.4 


0.7 


BP 


6 


22 


-94 


1.5-2.8 



TABLE I: Measured and calculated zero-bias differential con- 
ductances (in units of lO^^Co) for the different molecule-gold 
junctions under study. The computed spread reflects the low- 
est and highest conductance of the three geometries consid- 
ered for each junction (see SI for details). Experimental data 
for BP comes from Ref. J^. All other experimental data is 
taken from Ref. [22]. 



benzene diamine derivatives on flat gold (111) [36 . 

In |lj we report computed zero-bias conductances for 
DAS, APA, HDA, and BP junctions for both DFT-PBE 
and DFT+S, and compare with experiments [151 E2]- 
As with previous work [12||ig[4lH44], we find DFT-PBE 
overestimates measured low-bias conductances for DAS, 
APA, and BP junctions by more than an order of magni- 
tude. An exception is the HDA junction, where the over- 
estimate is just a factor of 4, consistent with Ref. [551 (56]. 
In contrast, DFT+S improves agreement with experi- 
ment significantly, predicting low-bias conductances to 
within a factor of two and shifting frontier orbital reso- 
nances to higher energies of 1.5 eV (BP-LUMO) or more 
away from Ep. For DAS, APA, and BP junctions, the to- 
tal correction Sqp + Scorr opens the gap between junc- 
tion HOMO and LUMO molecular resonances by more 
than 2.5 eV. For HDA, we find that Sqp and Scorr 
are of the same magnitude but opposite sign, resulting 
in a modest net correction of 0.15 eV, explaining its rel- 
atively accurate DFT-PBE zero-bias conductance by a 
cancellation of errors. 

Calculated IV characteristics for each junction are 
shown in [l] and compared with experiments. We first 
generate a DFT-PBE steady-state charge density, in- 
cluding a real-axis integration of the scattering-states in 
the bias window, at each bias voltage; the current is 
then determined from an integration of T{oj;V)^ sub- 
sequently computed on a dense energy grid, as described 
above. Then, in one approach, which we refer to as DFT- 
PBE(V), T{uj; V) is generated from DFT-PBE junction 
electronic structure, which is known to overestimate the 
zero-bias conductance; in the other approach, referred to 
as DFT(V)+S, T {oo; V) is generated from DFT+S junc- 
tion electronic structure, with S determined at zero-bias. 

As can be seen in [ij DFT(V)+S leads to excellent 
quantitative agreement with experiment, whereas the 
DFT-PBE (V) currents are too large by an order of mag- 
nitude (or more). Interestingly, integration of the zero- 
bias DFT-PBE transmission function (gray curves), as- 
suming a uniform potential drop across the junction, re- 
sults in a good estimate of the DFT-PBE(V) IV charac- 
teristics for DAS, APA, and HDA junctions, suggesting 



that T{uj;V) is weakly dependent on applied bias for 
these junctions. A similar integration results in an over- 
estimate by a factor of 8 to 10 for the BP junctions, which 
we explain below. 

In [2j we plot the evolution of the DFT-PBE charge 
density, Ap, and electrostatic Hartree potential, AUh^ 
with bias. Ap {V) and AUh {V) are both computed rel- 
ative to zero-bias densities and potentials, respectively. 
We compare this evolution to that generated assuming 
linear responses of the density and the Hartree poten- 



tial, i.e. Ap 



dv 



X V and AUh 



dV 



v=o 



respectively. We find that all 4 junctions remain in the 
linear response regime for the experimental range of bias 
(±0.8 V). This finding is consistent with the fact that 
the DFT+S molecular resonances are far from the bias 
window in these junctions, even for the highest biases 
achieved, leading to modest changes in occupation num- 
bers. Indeed, a simple estimate of the change in charge 
density under bias, obtained by fitting the DFT(V=0)+i; 
transmission function to a symmetric Lorentzian model 
for all four junctions (see SI), leads to a maximum ex- 
pected change in occupation by only 0.007 e~ . 

Interestingly for BP junctions ([2|, the DFT-PBE (V) 
densities and Hartree potentials remain in the linear re- 
sponse regime to 1 V, even though the DFT-PBE level 
alignment erroneously places the LUMO resonance in 
the bias window at this voltage. This indicates that, 
despite the incorrect DFT-PBE energy level alignment, 
self-consistent DFT-PBE (V) results are nonetheless in 
line with the DFT+S level alignment for the prediction 
of the densities under bias and do not show significant 
charging or any (non-physical) deviation from linear re- 
sponse. 

In [3j we report the evolution of the molecular reso- 
nances energies with bias for two BP junctions. For a 
symmetric junction, the resonance energies vary as F^, 
as expected from a quadratic Stark effect, with a maxi- 
mum shift (relative to the average chemical potential) of 
less than 0.06 eV for DAS and APA, and up to 0.2 eV 
for HDA and BP (|3] (a)). The small magnitude of these 
shifts explains why a simple integration of the zero-bias 
T {uj;V = 0) is so effective for these junctions. However, 
if the junctions are more asymmetric, such as the BP 
junction shown in[3](b), the resonance energies may vary 
more substantially. In fact, for the asymmetric BP junc- 
tion shown in|3](b), the LUMO resonance energy shifts 
linearly with V, with small second-order corrections to 
the energy. Appreciable first-order corrections to the 
wavefunction are also evident in this case, via a 50% re- 
duction in the resonant peak maximum. Integration of 
the zero-bias T(cj; V = 0) is clearly much less effective 
in this case, as T {uj; V) is varying significantly with V. 

Since the Stark shifts of the LUMO resonance for BP 
junctions are well-described by a simple function of V, 
a modified self-energy correction of the form S (V) = 
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FIG. 2: Junction structure, and DFT-PBE differences in electronic density and Hartree potentials in a BP junction for biases 
up to 1000 mV. (a) Junction structure. The red shaded areas indicate where densities and potentials are fixed to their bulk 
values; the frame in (a) indicates the region plotted in (b) and (c). (b) Ap and AUh, determined via self-consistent DFT- 
PBE(V) calculations, (c) Linear extrapolation result using p and Uh at zero bias and ^ and calculated at —100 mV. The 
color map represents differences in density, and the contour lines indicate isovalues of the change in Hartree potential scaled 
by the applied bias, i. e. = -50%, -40% • • • + 50%. 
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FIG. 3: Molecular resonance energies and IV character- 
istics for a (a) symmetric and (b) asymmetric BP junc- 
tion. Top: Calculated I{V) using self-consistent DFT(V)+S 
(black dots), DFT(V=0)+S (blue dotted hne), and one- 
shot DFT(V=0)+E(y) with both first-order (red crosses) 
and second-order (cyan squares) corrections to the zero bias 
Hamiltonian. Insets: Structure of the junctions. Bottom: 
Frontier orbital resonance energies with respect to the aver- 
age chemical potential of the junction as a function of applied 
voltage, self-consistently calculated with DFT(V)+S (dots), 
and best polynomial fits (line) of the change in resonance en- 
ergy with V, i.e. (a) AE{V) = -0.04^ + 0.2V^ and (b) 
AE{V) = -0.32^ + 0.02V^. 



I;gp + 5]corr + (V), with AE (V) from a polynomial 
fit of resonance energy shifts from just a few finite bias 
calculations, results in predicted currents in close agree- 
ment with those computed self-consistently at each V. 
Since our "Stark-corrected" self-energy ^{V) is calcula- 
ble at very low-bias or from the zero-bias polarizability, 



it can be used for efficient and quantitative prediction 
of nonlinear IVs and rectification ratios, at least for off- 
resonance tunneling through non-degenerate levels [61 . 

Our results for BP junctions, shown in Figures [2] and [3) 
are consistent with our findings for DAS, APA, and HDA. 
For each junction, we can draw the following conclusions. 
First, extrapolations using the zero-bias response func- 



tions 



dV 



and 



dV 



v=o 



result in accurate densities 



and potentials at other experimentally achievable biases, 
obviating the need in these cases to do laborious self- 
consistent calculations at many bias voltages. Second, 
self-consistency at DFT-PBE level is apparently suffi- 
cient for these junctions, despite errors in level align- 
ment, provided that the actual resonances are far from 
the bias window. Third, a straightforward Stark correc- 
tion AE (y) to the DFT-PBE resonance peak energy at 
zero-bias E{V = 0)^ along with a zero-bias self-energy 
correction, yields accurate IV characteristics, even for 
asymmetric junctions. 

The success of DFT+S at finite-bias is noteworthy and 
deserves further comment. Because of the small polariz- 
abilities of the molecules considered here, modifications 
of Sqp in relevant electric fields are just 0.01 eV for 
DAS and APA, 0.02 eV for HDA, and -0.002 eV for BP, 
a mere few % correction to the zero-bias self-term [62]. 
Moreover, using the DFT+S frontier orbital energies, the 
weak densities of states in the bias window for all junc- 
tions give rise to negligible changes in occupation under 
bias, smaller than 0.007 e~, implying equally negligible 
changes in Scorr for these junctions in the off-resonance 
limit. Thus corrections to S = 0) at finite bias will 
be extremely small. We expect the good agreement be- 
tween the DFT+S IV characteristics and experiments 
in [l] to deteriorate for more polarizable molecules (due 
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to the changes in the Hqp term), stronger couphng of 
resonances to the lead states (via a breakdown of the 
DFT+S approximation), and for molecular junctions be- 
yond the linear response regime of the density (associated 
with changes in Scorr term, and high-order terms in V 
in AE{V))^ e.g. for molecular resonances closer to the 
Fermi energy of the metal. 

To conclude, we have computed IV characteristics for 
four different molecular junctions, and compared them 
with STM-BJ experiments. While state-of-the-art meth- 
ods based on DFT-PBE largely overestimate measured 
currents, our pragmatic first-principles approach based 
on one-shot self-energy corrections that improve level 
alignment in the junction leads to excellent quantita- 
tive agreement. This method opens up new avenues for 
computationally inexpensive and quantitative modeling 
of non-equilibrium steady-state properties of molecular 
junctions: our results suggest that the finite-bias trans- 
mission function for both symmetric and asymmetric 
molecular junctions can be well- approximated by calcu- 
lating the changes in zero-bias resonance energies under 
small perturbative fields. These findings offer the possi- 
bility to more rapidly understand and predict functional 
nonlinear properties of junctions, such as the rectifica- 
tion, accurately and efficiently. 
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